EBEP1

library(readxl)
library(writexl)
library(tidyr)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ purrr     1.0.2
## ✔ forcats   1.0.0     ✔ readr     2.1.5
## ✔ ggplot2   3.5.1     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(dplyr)
library(ggplot2)
library(gtsummary)
library(ggforce)
library(stringr)
library(UpSetR)
library(forcats)
library(broom)

library(scales)  # pour percent_format()

## 
## Attachement du package : 'scales'
## 
## L'objet suivant est masqué depuis 'package:purrr':
## 
##     discard
## 
## L'objet suivant est masqué depuis 'package:readr':
## 
##     col_factor

library(ggtext)
library(GGally)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

library(UpSetR)
library(RColorBrewer)
library(purrr)


my_theme <- 
  theme( 
  panel.background = element_rect(fill = "white"),
  axis.line = element_line(color = "grey90"),
  legend.position = "top",
  legend.title = element_blank(),
  strip.background = element_rect(fill = "grey30"),
  strip.text = element_text(colour = "white"),
  text = element_text(face = "bold")
)

1. Description

Correspondance_EBEP1_2025 <-
  readxl::read_excel("./Raw_data/Correspondance_EBEP1_2025.xlsx")


Recode_of_every_possible_answer <-
  readxl::read_excel("./Raw_data/Recode_of_every_possible_answer.xlsx")

Reponse_EBE_P1_pre_filtered <-
  readxl::read_excel("./Raw_data/Survey_result.xlsx")

Reponse_EBE_P1 <- 
  Reponse_EBE_P1_pre_filtered[, (colnames(Reponse_EBE_P1_pre_filtered) %in% 
                                   Correspondance_EBEP1_2025$Nom)]

# Créer un vecteur de correspondance entre les noms et les codes
correspondance <- setNames(Correspondance_EBEP1_2025$Code, Correspondance_EBEP1_2025$Nom)

# Copier le dataframe d'origine
Reponse_EBE_P1_recode <- Reponse_EBE_P1

# Remplacer les noms de colonnes si correspondance trouvée
colnames(Reponse_EBE_P1_recode) <- ifelse(
  colnames(Reponse_EBE_P1_recode) %in% names(correspondance),
  correspondance[colnames(Reponse_EBE_P1_recode)],
  colnames(Reponse_EBE_P1_recode)  # laisser inchangé sinon
)


# Fonction de recodage d'une colonne
recode_2024 <- function(x) {
  x <- as.character(x)  # au cas où ce serait numérique
  x[x == "2024.0"] <- "2024"
  return(x)
}

# Appliquer à tout le dataframe (par exemple Reponse_EBE_P1)
Reponse_EBE_P1_recode <- as.data.frame(lapply(Reponse_EBE_P1_recode, recode_2024), stringsAsFactors = FALSE)

# # Step 1: Get all unique answers from Reponse_EBE_P1
# all_answers <- unique(unlist(Reponse_EBE_P1))
# 
# # Step 2: Check which answers are NOT in the recode list
# unknown_answers <- setdiff(all_answers, Recode_of_every_possible_answer$reponse_possible)
# 
# if(length(unknown_answers) > 0) {
#   cat("These answers have no recoding (not found in recode list):\n")
#   print(unknown_answers)
# } else {
#   cat("All answers found in the recode list. Proceeding with recoding...\n")
# }

# Step 3: Create a named vector for recoding for fast lookup
recode_vector <- 
  setNames(Recode_of_every_possible_answer$translation, 
           Recode_of_every_possible_answer$reponse_possible)

# Step 4: Replace answers by translation in the whole dataframe
# This will replace only known answers; unknown will become NA or stay same if you prefer

Reponse_EBE_P1_final <- as.data.frame(
  lapply(Reponse_EBE_P1_recode, function(col) {
    # Replace values by recode_vector if found, otherwise keep original or NA
    recoded_col <- recode_vector[as.character(col)]
    # Optional: keep original if not found
    recoded_col[is.na(recoded_col)] <- as.character(col)[is.na(recoded_col)]
    return(recoded_col)
  }),
  stringsAsFactors = FALSE
)
as_tibble(Reponse_EBE_P1_final)

## # A tibble: 5,443 × 221
##    University      Gender Age             Baccalaureate_Year Baccalaureate_Field
##    <chr>           <chr>  <chr>           <chr>              <chr>              
##  1 Montpellier     Female between 18 and… 2024.0             <NA>               
##  2 Antilles-Guyane Female between 18 and… 2024.0             <NA>               
##  3 Antilles-Guyane Male   between 18 and… 2024.0             <NA>               
##  4 Antilles-Guyane Female between 18 and… Between 2021 and … <NA>               
##  5 Antilles-Guyane Female between 18 and… 2024.0             <NA>               
##  6 Antilles-Guyane Female between 18 and… 2024.0             <NA>               
##  7 Antilles-Guyane Female between 18 and… 2024.0             <NA>               
##  8 Antilles-Guyane Female between 18 and… 2024.0             <NA>               
##  9 Antilles-Guyane Female between 18 and… 2024.0             <NA>               
## 10 Antilles-Guyane Female between 18 and… Between 2021 and … <NA>               
## # ℹ 5,433 more rows
## # ℹ 216 more variables: HighSchool_Last_Year_specialty_History_Geography <chr>,
## #   HighSchool_Last_Year_specialty__Humanities_Litterature_Philosophy <chr>,
## #   HighSchool_Last_Year_specialty_Humanities_Litterature_Languages <chr>,
## #   HighSchool_Last_Year_specialty_Maths <chr>,
## #   HighSchool_Last_Year_specialty_Physics <chr>,
## #   HighSchool_Last_Year_specialty_Earth_Sciences_and_Life <chr>, …

Reponse_EBE_P1_final <- as.data.frame(lapply(Reponse_EBE_P1_final, recode_2024), stringsAsFactors = FALSE)

Reponse_EBE_P1_final %>%
  writexl::write_xlsx("./Preprocessed_data/Reponse_EBE_P1_final.xlsx")

HAD_correspondance <-
  Correspondance_EBEP1_2025 %>%
  filter(Echelle == "HAD") 

HAD_correspondance$Nom_anglais <-
  str_remove_all(HAD_correspondance$Nom_anglais, "\\[|\\]")
HAD_correspondance$Nom_anglais <-
  str_remove_all(HAD_correspondance$Nom_anglais, "HAD ")
HAD_correspondance$Nom_anglais <-
  str_replace_all(HAD_correspondance$Nom_anglais, '"wound up"', 'wound up')


HAD_Reponse_EBE_P1_final <-
  Reponse_EBE_P1_final %>%
  select(contains("HAD"))

# Score Calcul ####
HAD_Reponse_EBE_P1_final_for_score <- HAD_Reponse_EBE_P1_final


# Plot HAD ####
# Récupérer les noms actuels des colonnes
old_names <- colnames(HAD_Reponse_EBE_P1_final)

# Créer un vecteur de nouveaux noms, en cherchant correspondance dans HAD_correspondance
new_names <- sapply(old_names, function(x) {
  idx <- which(HAD_correspondance$Code == x)
  if(length(idx) == 1) {
    return(HAD_correspondance$Nom_anglais[idx])
  } else {
    # Si pas de correspondance, on garde le nom original
    return(x)
  }
})

# Appliquer les nouveaux noms
colnames(HAD_Reponse_EBE_P1_final) <- new_names


HAD_Reponse_EBE_P1_final_long <- HAD_Reponse_EBE_P1_final %>%
    pivot_longer(
      cols = everything(),
      names_to = "Question",
      values_to = "Response"
    )
  
HAD_Reponse_EBE_P1_final_summary <- HAD_Reponse_EBE_P1_final_long %>%
    group_by(Question, Response) %>%
    summarise(n = n(), .groups = "drop") %>%
    group_by(Question) %>%
    mutate(percent = n / sum(n))
  
# Forcer l'ordre des modalités
HAD_Reponse_EBE_P1_final_summary <- HAD_Reponse_EBE_P1_final_summary %>%
  mutate(Response = factor(Response, levels = c("3 Most of the time", "2 Often", "1 From time to time", "0 Never")))

HAD_Reponse_EBE_P1_final_summary$Response %>% levels()

## [1] "3 Most of the time"  "2 Often"             "1 From time to time"
## [4] "0 Never"

# Calculer la somme de 0 et 1 par question pour mettre en gras si > 50%
bold_questions <- HAD_Reponse_EBE_P1_final_summary %>%
  filter(Response %in% c("0 Never", "1 From time to time")) %>%
  group_by(Question) %>%
  summarise(sum_percent = sum(percent)) %>%
  filter(sum_percent > 0.5) %>%
  pull(Question)

# Ajouter une colonne pour label en gras ou non
HAD_Reponse_EBE_P1_final_summary <- HAD_Reponse_EBE_P1_final_summary %>%
  mutate(bold_label = ifelse(Question %in% bold_questions, TRUE, FALSE))

# Préparer les labels pour afficher les % dans les barres (en % arrondi)
HAD_Reponse_EBE_P1_final_summary <- HAD_Reponse_EBE_P1_final_summary %>%
  mutate(percent_label = ifelse(percent > 0.02, paste0(round(percent*100), "%"), ""))  # n'affiche que si > 2%



# Define a gradient palette from light to dark blue (for example)
my_palette_HAD <- scales::gradient_n_pal(c("lightblue", "dodgerblue", "blue", "darkblue"))

# Generate 4 colors evenly spaced on the gradient
colors_HAD <- my_palette_HAD(seq(0, 1, length.out = 4))


ggplot(HAD_Reponse_EBE_P1_final_summary, aes(x = Question, y = percent, fill = Response)) +
  geom_bar(stat = "identity", position = "fill") +
  geom_text(aes(label = percent_label),
            position = position_fill(vjust = 0.5),
            color = "white",
            size = 3) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  coord_flip() +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = colors_HAD) +
  labs(title = "Distribution of HAD responses by question",
       x = "Question",
       y = "Percentage",
       fill = "Response") +
  theme_minimal() 

Due to a technical error, we cannot compute a score for HAD.

# Filtrer la correspondance pour l'échelle Warwick-Edinburgh
Warwick_correspondance <- Correspondance_EBEP1_2025 %>%
  filter(Echelle == "Warwick_Edinburgh")

# Sélectionner uniquement les colonnes liées à Warwick dans les réponses
WE_Reponse_EBE_P1_final <- Reponse_EBE_P1_final %>%
  select(contains("Warwick_Edinburgh"))


WE_Reponse_EBE_P1_final_for_score <- WE_Reponse_EBE_P1_final %>%
  mutate(across(everything(), ~ {
    # Retirer tous les caractères sauf chiffres, signes négatifs et points (virgules selon besoin)
    cleaned <- str_replace_all(as.character(.), "[^0-9.-]", "")
    # Convertir en numérique
    as.numeric(cleaned)
  }))


WE_Reponse_EBE_P1_final_score <-
  WE_Reponse_EBE_P1_final_for_score %>% rowSums() %>% as.data.frame()

colnames(WE_Reponse_EBE_P1_final_score) <- "Warwick_Edinburg_total"

# Renommage des colonnes avec noms anglais (si disponibles)
old_names <- colnames(WE_Reponse_EBE_P1_final)

new_names <- sapply(old_names, function(x) {
  idx <- which(Warwick_correspondance$Code == x)
  if(length(idx) == 1) {
    return(Warwick_correspondance$Nom_anglais[idx])
  } else {
    return(x)
  }
})

colnames(WE_Reponse_EBE_P1_final) <- new_names

# Passage en format long
WE_Reponse_EBE_P1_final_long <- WE_Reponse_EBE_P1_final %>%
  pivot_longer(
    cols = everything(),
    names_to = "Question",
    values_to = "Response"
  )

# Résumé des réponses
WE_Reponse_EBE_P1_final_summary <- WE_Reponse_EBE_P1_final_long %>%
  group_by(Question, Response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(Question) %>%
  mutate(percent = n / sum(n))


# Définir l’ordre des modalités de réponse typiques WEMWBS (à adapter selon ton codage réel)


WE_Reponse_EBE_P1_final_summary <- WE_Reponse_EBE_P1_final_summary %>%
  mutate(Response = factor(Response,
                           levels = c("5 All the time", "4 Often", "3 Sometimes", "2 Rarely", "1 never")))


WE_Reponse_EBE_P1_final_summary <- WE_Reponse_EBE_P1_final_summary %>%
  mutate(bold_label = ifelse(Question %in% bold_questions, TRUE, FALSE),
         percent_label = ifelse(percent > 0.02, paste0(round(percent*100), "%"), ""))

# Palette de couleurs pour WEMWBS
my_palette_WE <- scales::gradient_n_pal(c("lightgreen", "mediumseagreen", "green", "seagreen", "darkgreen"))
colors_WE <- my_palette_WE(seq(0, 1, length.out = 5))

# Visualisation
ggplot(WE_Reponse_EBE_P1_final_summary, aes(x = Question, y = percent, fill = Response)) +
  geom_bar(stat = "identity", position = "fill") +
  geom_text(aes(label = percent_label),
            position = position_fill(vjust = 0.5),
            color = "white",
            size = 3) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  coord_flip() +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = colors_WE) +
  labs(title = "Distribution of Warwick-Edinburgh responses by question",
       x = "Question",
       y = "Percentage",
       fill = "Response") +
  theme_minimal()

# 1. Filtrer la correspondance SCOFF
SCOFF_correspondance <- Correspondance_EBEP1_2025 %>%
  filter(Echelle == "SCOFF")

# 2. Sélectionner les colonnes SCOFF dans les réponses
SCOFF_Reponse_EBE_P1_final <- Reponse_EBE_P1_final %>%
  select(contains("SCOFF"))

# 3. Nettoyage et transformation en score binaire (0/1)
SCOFF_Reponse_EBE_P1_final_for_score <- SCOFF_Reponse_EBE_P1_final %>%
  mutate(across(everything(), ~ case_when(
    str_to_lower(.) %in% c("yes", "oui") ~ 1,
    str_to_lower(.) %in% c("no", "non") ~ 0,
    TRUE ~ NA_real_  # ou NA_integer_ si tu veux des entiers
  )))

# 4. Calcul du score total par individu
SCOFF_Reponse_EBE_P1_final_score <- SCOFF_Reponse_EBE_P1_final_for_score %>%
  rowSums(na.rm = TRUE) %>%
  as.data.frame()
colnames(SCOFF_Reponse_EBE_P1_final_score) <- "SCOFF_total"

# 5. Application du cutoff clinique
SCOFF_Reponse_EBE_P1_final_score <- SCOFF_Reponse_EBE_P1_final_score %>%
  mutate(SCOFF_risque_TCA = ifelse(SCOFF_total >= 2, "Risk_SCOFF", "No_Risk"))

# Barplot

# Préparation du résumé SCOFF
SCOFF_summary <- SCOFF_Reponse_EBE_P1_final_score %>%
  group_by(SCOFF_total) %>%
  summarise(n = n(), .groups = "drop") %>%
  mutate(percent = n / sum(n),
         SCOFF_total = factor(SCOFF_total, levels = sort(unique(SCOFF_total))),
         SCOFF_num = as.numeric(as.character(SCOFF_total)) # Conversion en numérique
  ) %>%
  mutate(SCOFF_total = factor(SCOFF_total, levels = as.character(0:5)))

# Calculer la hauteur cumulée au seuil SCOFF_num == 2
cutoff_y <- sum(SCOFF_summary$percent[SCOFF_summary$SCOFF_num < 2])

ggplot(SCOFF_summary, aes(x = 1, y = percent, fill = SCOFF_num)) +
  geom_bar(stat = "identity", width = 0.7, position = "stack") +
  geom_text(aes(label = paste0(round(percent * 100), "%")),
            position = position_stack(vjust = 0.5),
            color = "white", size = 5) +
  scale_fill_gradientn(
    colors = c("#FFC0CB", "#DA70D6", "#8A2BE2"),  # dégradé rose -> violet
    breaks = 0:5,
    labels = 0:5,
    limits = c(0,5),
    name = "SCOFF Total Score"
  ) +
  # Ligne horizontale à la hauteur du cutoff
  geom_hline(yintercept = cutoff_y, linetype = "dashed", color = "black", linewidth = 1) +
  
  # Annotation "Not at risk" à gauche sous la ligne
  annotate("text", x = 0.7, y = cutoff_y / 2, label = "Not at risk", angle = 0, size = 5, hjust = 1) +
  
  # Annotation "At risk" à gauche au-dessus de la ligne
  annotate("text", x = 0.7, y = cutoff_y + (1 - cutoff_y) / 2, label = "At risk", angle = 0, size = 5, hjust = 1) +
  
  coord_flip() +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    title = "Distribution of SCOFF Total Scores",
    x = "",
    y = "Proportion"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size = 12),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    legend.position = "bottom",
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 11)
  )

# 6. Renommage des colonnes avec noms anglais (si disponibles)
old_names <- colnames(SCOFF_Reponse_EBE_P1_final)
new_names <- sapply(old_names, function(x) {
  idx <- which(SCOFF_correspondance$Code == x)
  if(length(idx) == 1) {
    return(SCOFF_correspondance$Nom_anglais[idx])
  } else {
    return(x)
  }
})
colnames(SCOFF_Reponse_EBE_P1_final) <- new_names

# 7. Passage en format long
SCOFF_Reponse_EBE_P1_final_long <- SCOFF_Reponse_EBE_P1_final %>%
  pivot_longer(
    cols = everything(),
    names_to = "Question",
    values_to = "Response"
  )

# 8. Résumé des réponses
SCOFF_Reponse_EBE_P1_final_summary <- SCOFF_Reponse_EBE_P1_final_long %>%
  group_by(Question, Response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(Question) %>%
  mutate(percent = n / sum(n))

# 9. Forcer l’ordre (0 = non, 1 = oui)
SCOFF_Reponse_EBE_P1_final_summary <- SCOFF_Reponse_EBE_P1_final_summary %>%
  mutate(Response = factor(Response, levels = c("Yes", "No"), labels = c("Yes", "No")))

# 10. Préparation des étiquettes
SCOFF_Reponse_EBE_P1_final_summary <- SCOFF_Reponse_EBE_P1_final_summary %>%
  mutate(percent_label = ifelse(percent > 0.02, paste0(round(percent*100), "%"), ""))

# 11. Palette et graphe
colors_SCOFF <- c("darkred", "darkgreen")

ggplot(SCOFF_Reponse_EBE_P1_final_summary, aes(x = Question, y = percent, fill = Response)) +
  geom_bar(stat = "identity", position = "fill") +
  geom_text(aes(label = percent_label),
            position = position_fill(vjust = 0.5),
            color = "white",
            size = 3) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  coord_flip() +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = colors_SCOFF) +
  labs(title = "Distribution of SCOFF responses by question",
       x = "Question",
       y = "Percentage",
       fill = "Answer") +
  theme_minimal()

# Filtrer la correspondance pour MBI-SS
MBI_correspondance <- Correspondance_EBEP1_2025 %>%
  filter(Echelle == "French_MBI_SS")

# Sélectionner uniquement les colonnes MBI dans les réponses
MBI_Reponse_EBE_P1_final <- Reponse_EBE_P1_final %>%
  select(contains("French_MBI_SS"))

MBI_Reponse_EBE_P1_final_for_score <- MBI_Reponse_EBE_P1_final %>%
  mutate(across(everything(), ~ {
    cleaned <- str_replace_all(as.character(.), "[^0-9.-]", "")
    as.numeric(cleaned)
  }))


MBI_correspondance$Nom_anglais <- str_remove_all(MBI_correspondance$Nom_anglais, "\\[MBI-SS\\] ")


# #Epuismeent Emotionnel 
# MBI_Reponse_EBE_P1_final_for_score$French_MBI_SS_emotionally_drained
# MBI_Reponse_EBE_P1_final_for_score$French_MBI_SS_feeling_exhausted_day_end
# MBI_Reponse_EBE_P1_final_for_score$French_MBI_SS_feeling_exhausted_studies
# MBI_Reponse_EBE_P1_final_for_score$French_MBI_SS_tired_morning_university
# MBI_Reponse_EBE_P1_final_for_score$French_MBI_SS_studying_stressful
# 
# 
# #Cynisme
# MBI_Reponse_EBE_P1_final_for_score$French_MBI_SS_less_interest_studies
# MBI_Reponse_EBE_P1_final_for_score$French_MBI_SS_less_enthusiastic_studies
# MBI_Reponse_EBE_P1_final_for_score$French_MBI_SS_more_cynical_about_studies
# MBI_Reponse_EBE_P1_final_for_score$French_MBI_SS_doubt_study_meaning
# 
# #Efficacité académique
# MBI_Reponse_EBE_P1_final_for_score$French_MBI_SS_problem_solving_effective
# MBI_Reponse_EBE_P1_final_for_score$French_MBI_SS_effective_contribution
# MBI_Reponse_EBE_P1_final_for_score$French_MBI_SS_consider_good_student
# MBI_Reponse_EBE_P1_final_for_score$French_MBI_SS_learned_interesting_things
# MBI_Reponse_EBE_P1_final_for_score$French_MBI_SS_feeling_energized_goals
# MBI_Reponse_EBE_P1_final_for_score$French_MBI_SS_feeling_effective_in_class


MBI_Reponse_EBE_P1_final_for_score <- MBI_Reponse_EBE_P1_final_for_score %>%
  mutate(
    # Épuisement émotionnel : somme des 5 items
    Emotional_Exhaustion = rowSums(select(., 
                                          French_MBI_SS_emotionally_drained,
                                          French_MBI_SS_feeling_exhausted_day_end,
                                          French_MBI_SS_feeling_exhausted_studies,
                                          French_MBI_SS_tired_morning_university,
                                          French_MBI_SS_studying_stressful), na.rm = TRUE),
    
    # Cynisme : somme des 4 items
    Cynicism = rowSums(select(., 
                              French_MBI_SS_less_interest_studies,
                              French_MBI_SS_less_enthusiastic_studies,
                              French_MBI_SS_more_cynical_about_studies,
                              French_MBI_SS_doubt_study_meaning), na.rm = TRUE),
    
    # Efficacité académique : somme des 6 items
    Academic_Efficacy = rowSums(select(., 
                                       French_MBI_SS_problem_solving_effective,
                                       French_MBI_SS_effective_contribution,
                                       French_MBI_SS_consider_good_student,
                                       French_MBI_SS_learned_interesting_things,
                                       French_MBI_SS_feeling_energized_goals,
                                       French_MBI_SS_feeling_effective_in_class), na.rm = TRUE),
    
    
    # Classification de la sévérité
    Emotional_Exhaustion_Level = case_when(
      Emotional_Exhaustion >= 5  & Emotional_Exhaustion <= 13 ~ "Low",
      Emotional_Exhaustion >= 14 & Emotional_Exhaustion <= 22 ~ "Moderate",
      Emotional_Exhaustion >= 23 & Emotional_Exhaustion <= 30 ~ "High",
      TRUE ~ NA_character_
    ),
    
    Cynicism_Level = case_when(
      Cynicism >= 4  & Cynicism <= 10 ~ "Low",
      Cynicism >= 11 & Cynicism <= 17 ~ "Moderate",
      Cynicism >= 18 & Cynicism <= 24 ~ "High",
      TRUE ~ NA_character_
    ),
    
    Academic_Efficacy_Level = case_when(
      Academic_Efficacy >= 6  & Academic_Efficacy <= 16 ~ "Low",
      Academic_Efficacy >= 17 & Academic_Efficacy <= 26 ~ "Moderate",
      Academic_Efficacy >= 27 & Academic_Efficacy <= 36 ~ "High",
      TRUE ~ NA_character_
    )
  )
MBI_Reponse_EBE_P1_final_score <-
  MBI_Reponse_EBE_P1_final_for_score %>%
  select(
    # Scores globaux en premier
    Emotional_Exhaustion,
    Emotional_Exhaustion_Level,
    
    
    # Items Épuisement Émotionnel
    French_MBI_SS_emotionally_drained,
    French_MBI_SS_feeling_exhausted_day_end,
    French_MBI_SS_feeling_exhausted_studies,
    French_MBI_SS_tired_morning_university,
    French_MBI_SS_studying_stressful,
    
    Cynicism,
    Cynicism_Level,
    # Items Cynisme
    French_MBI_SS_less_interest_studies,
    French_MBI_SS_less_enthusiastic_studies,
    French_MBI_SS_more_cynical_about_studies,
    French_MBI_SS_doubt_study_meaning,
    
    Academic_Efficacy,
    Academic_Efficacy_Level,
    # Items Efficacité Académique
    French_MBI_SS_problem_solving_effective,
    French_MBI_SS_effective_contribution,
    French_MBI_SS_consider_good_student,
    French_MBI_SS_learned_interesting_things,
    French_MBI_SS_feeling_energized_goals,
    French_MBI_SS_feeling_effective_in_class,
    
    # Niveaux de sévérité
    )

# Plot Sévérité


# Long format + regroupement pour les proportions
MBI_levels_long <- MBI_Reponse_EBE_P1_final_score %>%
  select(contains("Level")) %>%
  pivot_longer(cols = everything(),
               names_to = "Dimension",
               values_to = "Severity_Level") %>%
  mutate(
    Dimension = case_when(
      Dimension == "Emotional_Exhaustion_Level" ~ "Emotional Exhaustion",
      Dimension == "Cynicism_Level" ~ "Cynicism",
      Dimension == "Academic_Efficacy_Level" ~ "Academic Efficacy"
    )
  ) %>%
  count(Dimension, Severity_Level) %>%
  group_by(Dimension) %>%
  mutate(
    percent = n / sum(n),
    Severity_Level = factor(Severity_Level, levels = c("Low", "Moderate", "High")),
    Dimension = factor(Dimension, levels = c("Emotional Exhaustion", "Cynicism", "Academic Efficacy"))
  )
MBI_levels_long$Severity_Level <- factor(
  MBI_levels_long$Severity_Level,
  levels = c("High", "Moderate", "Low")
)


ggplot(MBI_levels_long, aes(x = Dimension, y = percent, fill = Severity_Level)) +
  geom_bar(stat = "identity", position = "stack", width = 0.7) +
  geom_text(aes(label = paste0(round(percent * 100), "%")),
            position = position_stack(vjust = 0.5),
            color = "white", size = 5) +
  coord_flip() +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(
    values = c(
      "Low" = "#91cf60",       # green
      "Moderate" = "#fdae61",  # orange
      "High" = "#d73027"     # red
    ),
    name = "Level"
  ) +
  labs(
    title = "Burnout Severity Levels by Dimension",
    x = "",
    y = "Proportion"
  ) +
  my_theme

# Plot de tout  ####

MBI_Reponse_EBE_P1_final_score_filtered <-
  MBI_Reponse_EBE_P1_final %>%
  dplyr::select(all_of(MBI_correspondance$Code))

new_colnames <- 
  MBI_correspondance$Nom_anglais[match(
    colnames(MBI_Reponse_EBE_P1_final_score_filtered), 
    MBI_correspondance$Code)]

colnames(MBI_Reponse_EBE_P1_final_score_filtered) <- new_colnames



# Step 1: Convert data to long format
MBI_long <- MBI_Reponse_EBE_P1_final_score_filtered %>%
  pivot_longer(
    cols = everything(),
    names_to = "Question",
    values_to = "Response"
  )

# Step 2: Summarize responses per question
MBI_summary <- MBI_long %>%
  group_by(Question, Response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(Question) %>%
  mutate(percent = n / sum(n))

# Step 3: Define the factor levels of responses in French validated MBI-SS
# Adapt these levels according to your real coding (example given below)
response_levels <- c(
  "6 Always",
  "5.0", 
  "4.0", 
  "3.0", 
  "2.0", 
  "1 Never"
)

MBI_summary <- MBI_summary %>%
  mutate(Response = factor(Response, levels = response_levels))

# Optional: add labels if needed, for example, highlight specific questions
# (You can define a vector `bold_questions` with question names to highlight)
bold_questions <- c() # fill with your important question names if any

MBI_summary <- MBI_summary %>%
  mutate(bold_label = ifelse(Question %in% bold_questions, TRUE, FALSE),
         percent_label = ifelse(percent > 0.02, paste0(round(percent*100), "%"), ""))

# Step 4: Define a color palette for your responses (7 levels)
my_palette_MBI <- scales::gradient_n_pal(c(
  "lightblue", "skyblue", "deepskyblue", 
  "dodgerblue", "blue", "darkblue", "navy"
))
colors_MBI <- my_palette_MBI(seq(0, 1, length.out = length(response_levels)))

# Step 5: Plot
ggplot(MBI_summary, aes(x = Question, y = percent, fill = Response)) +
  geom_bar(stat = "identity", position = "fill") +
  geom_text(aes(label = percent_label),
            position = position_fill(vjust = 0.5),
            color = "white",
            size = 3) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  coord_flip() +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = colors_MBI) +
  labs(title = "Distribution of French MBI-SS responses by question",
       x = "Question",
       y = "Percentage",
       fill = "Response") +
  theme_minimal()

All_scores <-
  cbind((MBI_Reponse_EBE_P1_final_score %>%
          dplyr::select(-contains("French"))), 
       (SCOFF_Reponse_EBE_P1_final_score),
       WE_Reponse_EBE_P1_final_score)


All_scores %>% dplyr::select(where(is.numeric)) %>% ggpairs()

# Count responses per university
df_univ <- Reponse_EBE_P1_final %>%
  count(University)


# Optional: sort by descending count
df_univ <- df_univ %>%
  arrange(desc(n))

# Plot
ggplot(df_univ, aes(x = reorder(University, n), y = n)) +
  geom_bar(stat = "identity", fill = "darkblue") +
  geom_text(aes(label = n),
            hjust = 1.1,
            color = "white",
            size = 3.5) +
  coord_flip() +
  scale_y_continuous(expand = c(0.01, 0)) +  # small space at the start
  labs(title = "Distribution of responses by University",
       x = "University",
       y = "Number of responses") +
  my_theme

# Préparer les données
Score_and_univeristy <- Reponse_EBE_P1_final %>%
  dplyr::select("University") %>%
  cbind(All_scores) %>%
  select(c("University", where(is.numeric)))

# Format long
data_long <- Score_and_univeristy %>%
  pivot_longer(
    cols = -University,
    names_to = "Score_Type",
    values_to = "Score"
  )

# Diviser les données par Score_Type
data_split <- data_long %>% group_split(Score_Type)
score_names <- data_long %>% distinct(Score_Type) %>% pull()

# Créer les plots
plots <- map2(data_split, score_names, function(df, score_name) {
  
  # Médiane nationale
  national_median <- df %>%
    summarize(National_Median = median(Score, na.rm = TRUE)) %>%
    pull(National_Median)
  
  # Statistiques par université : médiane, max, effectif
  university_stats <- df %>%
    group_by(University) %>%
    summarize(
      Uni_Median = median(Score, na.rm = TRUE),
      Max_Score = max(Score, na.rm = TRUE),
      N = n(),
      .groups = "drop"
    )
  
  # Première université pour placer le label national median
  first_uni <- df %>% pull(University) %>% unique() %>% sort() %>% first()
  
  ggplot(df, aes(x = University, y = Score, fill = University)) +
    geom_violin(alpha = 0.3, colour = "grey40", draw_quantiles = c(0.25, 0.5, 0.75)) +
    geom_sina(alpha = 0.3) +
    
    # Ligne médiane nationale
    geom_hline(yintercept = national_median, linetype = "dashed", color = "black") +
    
    # Texte médiane nationale à gauche (sur la première université)
    annotate(
      "text",
      x = first_uni,
      y = national_median,
      label = paste0("National median: ", round(national_median, 1)),
      hjust = 0, vjust = -1, color = "red3", size = 3
    ) +
    
    # Texte médiane + effectif placé juste au-dessus du max de chaque université
    geom_text(
      data = university_stats, 
      aes(x = University, y = Max_Score + 0.5, label = paste0("(n = ", N, ")")),
      vjust = 0, size = 3, color = "black", inherit.aes = FALSE
    ) +
    
    labs(title = paste("Score:", score_name),
         x = "University", y = "Score") +
    theme_minimal() +
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1),
      legend.position = "none"
    )
})
plots

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

make_upset_plot <- function(data, prefix, keep_order = FALSE) {
  # 1. Sélectionner et recoder en binaire
  df_binary <- data %>%
    select(contains(prefix)) %>%
    mutate(across(everything(), ~ ifelse(.x == "Yes", 1,
                                         ifelse(.x == "No", 0, NA)))) %>%
    filter(if_all(everything(), ~ !is.na(.)))
  
  # 2. Renommer les colonnes en retirant le préfixe
  colnames(df_binary) <- str_remove(colnames(df_binary), prefix)
  
  # 3. Calculer le nombre de "Yes" et trier les colonnes
  yes_counts <- colSums(df_binary)
  ordered_columns <- names(sort(yes_counts, decreasing = TRUE))
  df_binary <- df_binary[, ordered_columns]
  
  # 4. Ajouter "(n = X)" aux noms de colonnes
  new_colnames <- paste0(ordered_columns, " (n = ", yes_counts[ordered_columns], ")")
  colnames(df_binary) <- new_colnames
  
  # 5. Afficher le plot UpSet
  upset(
    df_binary,
    sets = new_colnames,
    order.by = "freq",
    keep.order = keep_order
  )
}

make_upset_plot(Reponse_EBE_P1_final, "HighSchool_Last_Year_specialty_", keep_order = TRUE)

make_upset_plot(Reponse_EBE_P1_final, "HighSchool_SecondLast_Year_specialty_", keep_order = TRUE)

bacc_single <- Reponse_EBE_P1_final %>%
  filter(!is.na(Baccalaureate_Field)) %>%
  count(Baccalaureate_Field) %>%
  mutate(freq = n / sum(n) * 100)

my_colors <- brewer.pal(4, "Set2")  # soft pastel colors, great for 4 categories

ggplot(bacc_single, aes(x = 1, y = freq, fill = Baccalaureate_Field)) +
  geom_bar(stat = "identity", width = 0.5) +
  geom_text(aes(label = paste0(round(freq, 1), "%")),
            position = position_stack(vjust = 0.5), size = 4, color = "black") +
  coord_flip() +
  scale_fill_manual(values = my_colors) +
  labs(
    title = "Distribution of Baccalaureate Fields (in %)",
    x = NULL,
    #y = "Percentage (%)",
    fill = "Field"
  ) +
  theme_minimal() +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())

make_upset_plot(Reponse_EBE_P1_final, "Info_source_", keep_order = TRUE)

make_upset_plot(Reponse_EBE_P1_final %>% dplyr::select(-Psychiatric_history_Yes), "Psychiatric_history_", keep_order = TRUE)

make_upset_plot(Reponse_EBE_P1_final, "Violence_suffered_from_", keep_order = TRUE)

make_upset_plot(Reponse_EBE_P1_final, "Violence_type_", keep_order = TRUE) 

make_upset_plot(Reponse_EBE_P1_final, "Tutoring_service_looked_for_", keep_order = TRUE) 

make_upset_plot(Reponse_EBE_P1_final, "Private_preparation_service_looked_for_", keep_order = TRUE) 

make_upset_plot(Reponse_EBE_P1_final, "Academic_pre_entry_", keep_order = TRUE) 

make_upset_plot(Reponse_EBE_P1_final %>%
                  dplyr::select(-Workplace_suitable_for_your_needs), "Workplace_", keep_order = TRUE) 

# make_upset_plot(Reponse_EBE_P1_final, "Suicidal_thoughts_Frequency_", keep_order = TRUE) 
#  Not used anymore

Reponse_EBE_P1_final$Fifth_choice_health_track <- Reponse_EBE_P1_final$Fifth_choice_health_track %>%
  dplyr::recode(
    "No fifth choice other than fourth" = "No fifth choice different from fourth"
  )



health_track_long <- Reponse_EBE_P1_final %>%
  select(contains("health_track")) %>%
  mutate(Respondent_ID = row_number()) %>%
  pivot_longer(
    cols = -Respondent_ID,
    names_to = "Choice_Order",
    values_to = "Health_Track"
  )



health_track_long <- Reponse_EBE_P1_final %>%
  select(contains("health_track")) %>%
  mutate(Respondent_ID = row_number()) %>%
  pivot_longer(
    cols = -Respondent_ID,
    names_to = "Choice_Order",
    values_to = "Health_Track"
  ) %>%
  mutate(
    # 3. Ordonner les niveaux de Choice_Order
    Choice_Order = factor(Choice_Order, levels = c(
      "First_choice_health_track",
      "Second_choice_health_track",
      "Third_choice_health_track",
      "Fourth_choice_health_track",
      "Fifth_choice_health_track"
    ))
  )

# 3. Plot
ggplot(health_track_long %>% filter(!is.na(Health_Track)),
       aes(x = fct_infreq(Health_Track))) +
  geom_bar(fill = "#3182bd") +
  geom_text(stat = "count", aes(label = after_stat(count)), hjust = -0.1, size = 3.5) +
  coord_flip() +
  facet_wrap(~ Choice_Order, scales = "free_y") +
  labs(
    title = "Distribution of Health Track Choices by Rank",
    x = "Health Track",
    y = "Number of Respondents"
  ) +
  my_theme +
  theme(strip.text = element_text(size = 12),
        plot.title = element_text(face = "bold")) +
  expand_limits(y = 1.05 * max(table(health_track_long$Health_Track, health_track_long$Choice_Order)))

# Step 1: Reshape to long format
  health_track_long_2 <- Reponse_EBE_P1_final %>%
  select(contains("health_track")) %>%
  mutate(Respondent_ID = row_number()) %>%
  pivot_longer(
    cols = -Respondent_ID,
    names_to = "Choice_Order",
    values_to = "Health_Track"
  ) %>%
  mutate(
    Health_Track = if_else(str_detect(Health_Track, "No .* choice different from .*"),
                           NA_character_,
                           Health_Track),
    Choice_Order = factor(Choice_Order,
                          levels = c("First_choice_health_track",
                                     "Second_choice_health_track",
                                     "Third_choice_health_track",
                                     "Fourth_choice_health_track",
                                     "Fifth_choice_health_track"),
                          labels = c("1st", "2nd", "3rd", "4th", "5th"))
  ) %>%
  filter(!is.na(Health_Track))

# Count number of respondents per Health_Track and Choice_Order
data_count <- health_track_long_2 %>%
  group_by(Health_Track, Choice_Order) %>%
  summarise(count = n(), .groups = "drop")

# Reorder health tracks by overall frequency
health_levels <- health_track_long_2 %>%
  count(Health_Track) %>%
  arrange(desc(n)) %>%
  pull(Health_Track)

data_count <- data_count %>%
  mutate(Health_Track = factor(Health_Track, levels = health_levels),
         Choice_Order = fct_rev(Choice_Order))  # Reverse choice order (5th -> 1st)

# Plot with counts inside bars
ggplot(data_count, aes(x = Health_Track, y = count, fill = Choice_Order)) +
  geom_col() +
  geom_text(aes(label = count), position = position_stack(vjust = 0.5), size = 3) +
  coord_flip() +
  scale_fill_brewer(palette = "Blues", direction = -1) +
  labs(
    title = "Number of Respondents by Health Track and Choice Order",
    x = "Health Track",
    y = "Number of Respondents",
    fill = "Choice Order"
  ) +
  theme_minimal(base_size = 14)

Reponse_EBE_P1_final %>%
  dplyr::select(
    Current_track,
    PASS_before_LAS,
    Current_study_year,
    LAS_chosen_first,
    LAS_major,
    LAS_major_first_choice,
    LAS_continue_if_failed,
    PASS_chosen_first,
    PASS_minor,
    PASS_minor_first_choice,
    PASS_continue_if_failed
  ) %>% tbl_summary()

Characteristic	N = 5,4431
Current_track
LAS, Health Access License	1,650 (30%)
LSPS, Licence Science for Health, unique model in Strasbourg, Reims or Lille Catholique	418 (7.7%)
PASS, Specific Health Access Pathway	3,375 (62%)
PASS_before_LAS	478 (29%)
Unknown	3,793
Current_study_year
LAS 1	946 (57%)
LAS 2	594 (36%)
LAS 3	110 (6.7%)
Unknown	3,793
LAS_chosen_first
First choice	771 (82%)
I initially wanted to join the PASS stream	175 (18%)
Unknown	4,497
LAS_major
Biology	471 (29%)
Health sciences/STAPS	302 (18%)
Humanities_SocialSciences_Law_and_Arts	517 (31%)
Mathematics	40 (2.4%)
Mechanics/Engineering sciences/Computer science/Electronics	19 (1.2%)
Other	170 (10%)
Physics/Physical Chemistry/Chemistry	131 (7.9%)
Unknown	3,793
LAS_major_first_choice
no	351 (21%)
Yes	1,299 (79%)
Unknown	3,793
LAS_continue_if_failed
1 Strongly disagree	423 (26%)
2 Somewhat disagree	256 (16%)
3 Neutral	265 (16%)
4 Somewhat agree	561 (34%)
5 Totally agree	145 (8.8%)
Unknown	3,793
PASS_chosen_first
First choice	3,327 (99%)
I originally wanted to join the LAS stream	48 (1.4%)
Unknown	2,068
PASS_minor
Biology	523 (15%)
Health sciences/STAPS	193 (5.7%)
Humanities_SocialSciences_Law_and_Arts	763 (23%)
Mathematics	109 (3.2%)
Mechanics/Engineering sciences/Computer science/Electronics	55 (1.6%)
Other	957 (28%)
Physics/Physical Chemistry/Chemistry	619 (18%)
STAPS	156 (4.6%)
Unknown	2,068
PASS_minor_first_choice	2,644 (78%)
Unknown	2,068
PASS_continue_if_failed
1 Strongly disagree	634 (19%)
2 Somewhat disagree	536 (16%)
3 Neutral	628 (19%)
4 Somewhat agree	1,363 (40%)
5 Totally agree	214 (6.3%)
Unknown	2,068
1 n (%)

Reponse_EBE_P1_for_summary_withoutscore <-
  Reponse_EBE_P1_final %>%
  dplyr::select(
    -contains("Psychoactive"),
    -contains("Psychotropic"),
    # -all_of(HAD_correspondance$Code), #Finally we keep it just not a score
    -contains("HighSchool_Last_Year_specialty_"), 
    -contains("HighSchool_SecondLast_Year_specialty_"), 
    -contains("Baccalaureate_Field"), 
    -contains("University"),
    -contains("Info_source_"),
    -contains("Violence_suffered_from_"),
   # -contains("Violence_type_"),
    -contains("Tutoring_service_looked_for_"),
    -contains("Private_preparation_service_looked_for_"),
   -matches("Academic_pre_entry_(?!.*No)", perl = TRUE),
    -contains("Workplace_"),
    -contains("Suicidal_thoughts_Frequency_"),
   -matches("health_track(?!.*First)", perl = TRUE),
  -matches("Psychiatric_history_(?!.*(Yes|Eating|Anxiety|Depression))", perl = TRUE),
  -matches("Violence_type_(?!.*(Moral|Discrimination))", perl = TRUE),
   -PASS_before_LAS,
    -Current_study_year,
    -LAS_chosen_first,
    -LAS_major,
    -LAS_major_first_choice,
    -LAS_continue_if_failed,
    -PASS_chosen_first,
    -PASS_minor,
    -PASS_minor_first_choice,
    -PASS_continue_if_failed
  ) %>%
  as_tibble() 

Reponse_EBE_P1_for_summary <-
  cbind(Reponse_EBE_P1_for_summary_withoutscore, All_scores)

Reponse_EBE_P1_for_summary$Age <- factor(
  Reponse_EBE_P1_for_summary$Age,
  levels = c("between 18 and 20 years", "between 21 and 23", "24 years and over")
)

Reponse_EBE_P1_for_summary$Psychiatric_history_Anxiety <- Reponse_EBE_P1_for_summary$Psychiatric_history_Anxiety %>%
  na_if("NA") %>%        # transforme "NA" (chaine) en NA (valeur manquante)
  as.factor()

Reponse_EBE_P1_for_summary$Psychiatric_history_Depression <- Reponse_EBE_P1_for_summary$Psychiatric_history_Depression %>%
  na_if("NA") %>%        # transforme "NA" (chaine) en NA (valeur manquante)
  as.factor()



Reponse_EBE_P1_for_summary$Psychiatric_history_Eating_disorder <- Reponse_EBE_P1_for_summary$Psychiatric_history_Eating_disorder %>%
  na_if("NA") %>%        # transforme "NA" (chaine) en NA (valeur manquante)
  as.factor()



Reponse_EBE_P1_for_summary$Psychiatric_history_Yes <- Reponse_EBE_P1_for_summary$Psychiatric_history_Yes %>%
  na_if("NA") %>%        # transforme "NA" (chaine) en NA (valeur manquante)
  as.factor()

Reponse_EBE_P1_for_summary$Distance_home_univversity <-
  Reponse_EBE_P1_for_summary$Distance_home_univversity %>% as.factor()

Reponse_EBE_P1_for_summary$Baccalaureate_Year <- factor(
  Reponse_EBE_P1_for_summary$Baccalaureate_Year,
  levels = c(
    "Before 2019 (included)",
    "2020 Covid year",
    "Between 2021 and 2023", "2024"
  )
)

Reponse_EBE_P1_for_summary$Baccalaureate_honors <- factor(
  Reponse_EBE_P1_for_summary$Baccalaureate_honors,
  levels = c(
    "No honours",
    "Quite good",
    "Good",
    "Very good",
    "Highest_honours_more_than_18out20"
  )
)

Reponse_EBE_P1_for_summary$Weekly_study_hours <- factor(
  Reponse_EBE_P1_for_summary$Weekly_study_hours,
  levels = c(
    "Less than 15 hours",
    "15 to 30 hours",
    "30 to 45 hours",
    "45 to 60 hours",
    "60 to 75 hours",
    "More than 75 hours"
  )
)

Reponse_EBE_P1_for_summary$Weekly_hours_leisure_time <- factor(
  Reponse_EBE_P1_for_summary$Weekly_hours_leisure_time,
  levels = c(
    "Less than 5 hours",
    "5 to 10 hours",
    "10 to 15 hours",
    "15 to 20 hours",
    "More than 20 hours"
  )
)

Reponse_EBE_P1_for_summary$All_nighter_study <- factor(
  Reponse_EBE_P1_for_summary$All_nighter_study,
  levels = c(
    "Never",
    "Daily",
    "Weekly",
    "Monthly"
  )
)

Reponse_EBE_P1_for_summary$Weight_change_since_academic_year_beginning <- factor(
  Reponse_EBE_P1_for_summary$Weight_change_since_academic_year_beginning,
  levels = c(
    "No difference",
    "a loss of between 5 and 10% of your body weight",
    "a loss of at least 10% of your body weight",
    "an increase of between 5 and 10% in your weight",
    "an increase of at least 10% in your weight"
  )
)

Reponse_EBE_P1_for_summary$Distance_home_univversity <- factor(
  Reponse_EBE_P1_for_summary$Distance_home_univversity,
  levels = c("less than 15 minutes",
             "between 15 and 30 minutes",
             "between 30 minutes and 1 hour",
             "between 1h and 1h30",
             "more than 1h30"),
  ordered = TRUE
)

Reponse_EBE_P1_for_summary$Marital_status <- factor(
  Reponse_EBE_P1_for_summary$Marital_status,
  levels = c("Single", "Couple", "Married or civil union", "Widowed")
)

Reponse_EBE_P1_for_summary <- Reponse_EBE_P1_for_summary %>%
  mutate(Tutoring_participation_influence = factor(
    Tutoring_participation_influence,
    levels = c("Very negative", "Negative", "Neutral", "Positive", "Very positive"),
    ordered = FALSE
  ))

Reponse_EBE_P1_for_summary <- Reponse_EBE_P1_for_summary %>%
  mutate(Private_preparation_influence = factor(
    Private_preparation_influence,
    levels = c("Very negative", "Negative", "Neutral", "Positive", "Very positive"),
    ordered = FALSE
  ))


Reponse_EBE_P1_summary <-
  Reponse_EBE_P1_for_summary %>% tbl_summary()

Reponse_EBE_P1_summary %>% as_tibble() %>%
  writexl::write_xlsx("./Result/Reponse_EBE_P1_summary.xlsx")


Reponse_EBE_P1_summary_gender <-
  Reponse_EBE_P1_for_summary %>% 
  filter(Gender != "Other") %>%
  tbl_summary(by = "Gender") %>%
  add_p() %>% add_q() %>% bold_p()

## There was an error in 'add_p()/add_difference()' for variable 'French_MBI_SS_learned_interesting_things', p-value omitted:
## Error in stats::fisher.test(c("5.0", "6 Always", "4.0", "5.0", "2.0", : FEXACT erreur 7(location). LDSTP=17640 est trop petit pour ce problème,
##   (pastp=409.795, ipn_0:=ipoin[itp=571]=863, stp[ipn_0]=410.404).
## Augmentez la taille de l’environnement de travail ou considérez l’utilisation de ‘simulate.p.value=TRUE’.

## There was an error in 'add_p()/add_difference()' for variable 'Tutoring_participation_influence', p-value omitted:
## Error in stats::fisher.test(structure(c(5L, 3L, 3L, NA, NA, NA, 4L, 5L, : FEXACT erreur 7(location). LDSTP=17760 est trop petit pour ce problème,
##   (pastp=379.49, ipn_0:=ipoin[itp=554]=261, stp[ipn_0]=379.609).
## Augmentez la taille de l’environnement de travail ou considérez l’utilisation de ‘simulate.p.value=TRUE’.

## add_q: Adjusting p-values with
## `stats::p.adjust(x$table_body$p.value, method = "fdr")`

Reponse_EBE_P1_summary_gender %>% as_tibble() %>%
  writexl::write_xlsx("./Result/Reponse_EBE_P1_summary_gender.xlsx")


Reponse_EBE_P1_summary_track <-
  Reponse_EBE_P1_for_summary %>% tbl_summary(by = "Current_track") %>%
  add_p() %>% add_q() %>% bold_p()

## There was an error in 'add_p()/add_difference()' for variable 'Gender', p-value omitted:
## Error in stats::fisher.test(c("Female", "Female", "Male", "Female", "Female", : FEXACT erreur 6.  LDKEY=588 est trop petit pour ce problème,
##   (ii := key2[itp=1048] = 72414417, ldstp=17640)
## Essayez d’augmenter la taille de l’environnement de travail et peut-être aussi 'mult'
## There was an error in 'add_p()/add_difference()' for variable 'Baccalaureate_Year', p-value omitted:
## Error in stats::fisher.test(structure(c(4L, 4L, 4L, 3L, 4L, 4L, 4L, 4L, : FEXACT error 501.
## The hash table key cannot be computed because the largest key
## is larger than the largest representable int.
## The algorithm cannot proceed.
## Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.
## There was an error in 'add_p()/add_difference()' for variable 'Weekly_hours_leisure_time', p-value omitted:
## Error in stats::fisher.test(structure(c(1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, : FEXACT error 501.
## The hash table key cannot be computed because the largest key
## is larger than the largest representable int.
## The algorithm cannot proceed.
## Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.
## There was an error in 'add_p()/add_difference()' for variable 'Study_impact_on_family_life', p-value omitted:
## Error in stats::fisher.test(c("2 Negative impact", "2 Negative impact", : FEXACT error 501.
## The hash table key cannot be computed because the largest key
## is larger than the largest representable int.
## The algorithm cannot proceed.
## Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.
## There was an error in 'add_p()/add_difference()' for variable 'French_MBI_SS_learned_interesting_things', p-value omitted:
## Error in stats::fisher.test(c("5.0", "6 Always", "4.0", "5.0", "2.0", : FEXACT error 501.
## The hash table key cannot be computed because the largest key
## is larger than the largest representable int.
## The algorithm cannot proceed.
## Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.
## There was an error in 'add_p()/add_difference()' for variable 'French_MBI_SS_feeling_energized_goals', p-value omitted:
## Error in stats::fisher.test(c("6 Always", "4.0", "4.0", "6 Always", "5.0", : FEXACT error 501.
## The hash table key cannot be computed because the largest key
## is larger than the largest representable int.
## The algorithm cannot proceed.
## Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.
## There was an error in 'add_p()/add_difference()' for variable 'Marital_status', p-value omitted:
## Error in stats::fisher.test(structure(c(2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, : FEXACT error 501.
## The hash table key cannot be computed because the largest key
## is larger than the largest representable int.
## The algorithm cannot proceed.
## Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.
## There was an error in 'add_p()/add_difference()' for variable 'Dependent_children', p-value omitted:
## Error in stats::fisher.test(c("0.0", "0.0", "0.0", "0.0", "0.0", "0.0", : FEXACT error 501.
## The hash table key cannot be computed because the largest key
## is larger than the largest representable int.
## The algorithm cannot proceed.
## Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.
## There was an error in 'add_p()/add_difference()' for variable 'Housing', p-value omitted:
## Error in stats::fisher.test(c("In a group, sharing a flat, as a couple, etc...", : FEXACT error 501.
## The hash table key cannot be computed because the largest key
## is larger than the largest representable int.
## The algorithm cannot proceed.
## Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.
## There was an error in 'add_p()/add_difference()' for variable 'Tutoring_participation_influence', p-value omitted:
## Error in stats::fisher.test(structure(c(5L, 3L, 3L, NA, NA, NA, 4L, 5L, : FEXACT[f3xact()] error: hash key 3e+09 > INT_MAX, kyy=2984, it[i (= nco = 4)]= 0.
## Rather set 'simulate.p.value=TRUE'
## 
## There was an error in 'add_p()/add_difference()' for variable 'Tutoring_participation_can_be', p-value omitted:
## Error in stats::fisher.test(c("Essential", "Useful", "Useful", "Useful", : FEXACT error 501.
## The hash table key cannot be computed because the largest key
## is larger than the largest representable int.
## The algorithm cannot proceed.
## Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.
## There was an error in 'add_p()/add_difference()' for variable 'Private_preparation_participation', p-value omitted:
## Error in stats::fisher.test(c(NA, "I make partial use of the services offered", : FEXACT erreur 6.  LDKEY=609 est trop petit pour ce problème,
##   (ii := key2[itp=615] = 24342343, ldstp=18270)
## Essayez d’augmenter la taille de l’environnement de travail et peut-être aussi 'mult'
## There was an error in 'add_p()/add_difference()' for variable 'Private_preparation_influence', p-value omitted:
## Error in stats::fisher.test(structure(c(NA, 4L, NA, NA, 3L, 5L, 5L, 5L, : FEXACT erreur 6.  LDKEY=608 est trop petit pour ce problème,
##   (ii := key2[itp=567] = 88063243, ldstp=18240)
## Essayez d’augmenter la taille de l’environnement de travail et peut-être aussi 'mult'
## There was an error in 'add_p()/add_difference()' for variable 'Student_job_needed', p-value omitted:
## Error in stats::fisher.test(c("No", "No", "No", "No", "No", "No", "No", : FEXACT error 501.
## The hash table key cannot be computed because the largest key
## is larger than the largest representable int.
## The algorithm cannot proceed.
## Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.
## There was an error in 'add_p()/add_difference()' for variable 'Social_and_solidarity_grocery_store', p-value omitted:
## Error in stats::fisher.test(c("No", "No", "No", "Yes, but I don't go there", : FEXACT error 501.
## The hash table key cannot be computed because the largest key
## is larger than the largest representable int.
## The algorithm cannot proceed.
## Reduce the workspace, consider using 'simulate.p.value=TRUE' or another algorithm.
## add_q: Adjusting p-values with
## `stats::p.adjust(x$table_body$p.value, method = "fdr")`

Reponse_EBE_P1_summary_track %>% as_tibble() %>%
  writexl::write_xlsx("./Result/Reponse_EBE_P1_summary_track.xlsx")



Reponse_EBE_P1_summary_univr <-
  cbind(Reponse_EBE_P1_for_summary, (Reponse_EBE_P1_final %>% dplyr::select("University"))) %>% tbl_summary(by = "University") %>%
  add_overall()

Reponse_EBE_P1_summary_univr %>% as_tibble() %>%
  writexl::write_xlsx("./Result/Reponse_EBE_P1_summary_univr.xlsx")

library(ordinal)

## 
## Attachement du package : 'ordinal'

## L'objet suivant est masqué depuis 'package:dplyr':
## 
##     slice

Reponse_EBE_P1_for_regression <-
  Reponse_EBE_P1_for_summary %>%
  dplyr::select(-contains("Parent"),
         -Dependent_children,
         -Scolarship_level,
         -Tutoring_participation_can_be,
         -Private_preparation_participation_can_be,
         -Social_and_solidarity_grocery_store
         )
Reponse_EBE_P1_for_regression$SCOFF_total <- as.factor(Reponse_EBE_P1_for_regression$SCOFF_total)


translation_map <- setNames(Recode_of_every_possible_answer$translation_num,
                            Recode_of_every_possible_answer$translation)

# Fonction de recodage qui garde les valeurs non trouvées
safe_recode <- function(x, map) {
  recoded <- map[as.character(x)]
  ifelse(is.na(recoded), x, recoded)
}



# Application sur tout le dataframe
Reponse_EBE_P1_for_regression_recoded <- Reponse_EBE_P1_for_regression %>%
  mutate(across(everything(), ~ safe_recode(., translation_map)))
Reponse_EBE_P1_for_regression_recoded$SCOFF_total <-
  as.factor(Reponse_EBE_P1_for_regression_recoded$SCOFF_total)

Reponse_EBE_P1_for_regression_recoded$Age <- factor(
  Reponse_EBE_P1_for_regression_recoded$Age,
  levels = c("between 18 and 20 years", "between 21 and 23", "24 years and over")
)

Reponse_EBE_P1_for_regression_recoded$Baccalaureate_Year <- factor(
  Reponse_EBE_P1_for_regression_recoded$Baccalaureate_Year,
  levels = c(
    "Before 2019 (included)",
    "2020 Covid year",
    "Between 2021 and 2023", "2024"
  )
)



Reponse_EBE_P1_for_regression_recoded$Baccalaureate_honors <- factor(
  Reponse_EBE_P1_for_regression_recoded$Baccalaureate_honors,
  levels = c(
    "No honours",
    "Quite good",
    "Good",
    "Very good",
    "Highest_honours_more_than_18out20"
  )
)

Reponse_EBE_P1_for_regression_recoded$Weekly_study_hours <- factor(
  Reponse_EBE_P1_for_regression_recoded$Weekly_study_hours,
  levels = c(
    "Less than 15 hours",
    "15 to 30 hours",
    "30 to 45 hours",
    "45 to 60 hours",
    "60 to 75 hours",
    "More than 75 hours"
  )
)

Reponse_EBE_P1_for_regression_recoded$Weekly_hours_leisure_time <- factor(
  Reponse_EBE_P1_for_regression_recoded$Weekly_hours_leisure_time,
  levels = c(
    "Less than 5 hours",
    "5 to 10 hours",
    "10 to 15 hours",
    "15 to 20 hours",
    "More than 20 hours"
  )
)

Reponse_EBE_P1_for_regression_recoded$All_nighter_study <- factor(
  Reponse_EBE_P1_for_regression_recoded$All_nighter_study,
  levels = c(
    "Never",
    "Daily",
    "Weekly",
    "Monthly"
  )
)

Reponse_EBE_P1_for_regression_recoded$Weight_change_since_academic_year_beginning <- factor(
  Reponse_EBE_P1_for_regression_recoded$Weight_change_since_academic_year_beginning,
  levels = c(
    "No difference",
    "a loss of between 5 and 10% of your body weight",
    "a loss of at least 10% of your body weight",
    "an increase of between 5 and 10% in your weight",
    "an increase of at least 10% in your weight"
  )
)

Reponse_EBE_P1_for_regression_recoded$Distance_home_univversity <- factor(
  Reponse_EBE_P1_for_regression_recoded$Distance_home_univversity,
  levels = c("less than 15 minutes",
             "between 15 and 30 minutes",
             "between 30 minutes and 1 hour",
             "between 1h and 1h30",
             "more than 1h30"),
  ordered = TRUE
)


# Remplacer "Do not wish to answer" par NA
Reponse_EBE_P1_for_regression_recoded$Violence_suffered <- as.character(Reponse_EBE_P1_for_regression_recoded$Violence_suffered)
Reponse_EBE_P1_for_regression_recoded$Violence_suffered[Reponse_EBE_P1_for_regression_recoded$Violence_suffered == "Do not wish to answer"] <- NA

# Reconvertir en facteur (optionnel, selon ce que tu veux faire ensuite)
Reponse_EBE_P1_for_regression_recoded$Violence_suffered <- factor(
  Reponse_EBE_P1_for_regression_recoded$Violence_suffered,
  levels = c("No", "Yes")
)

# Find columns whose names contain "Violence_type"
violence_cols <- grep("Violence_type", names(Reponse_EBE_P1_for_regression_recoded), value = TRUE)

# Loop through these columns and replace "Do not wish to answer" by NA,
# then relevel factors with just "No" and "Yes" (dropping NA level)
for(col in violence_cols) {
  # Convert to character first
  Reponse_EBE_P1_for_regression_recoded[[col]] <- as.character(Reponse_EBE_P1_for_regression_recoded[[col]])
  
  # Replace "Do not wish to answer" by NA
  Reponse_EBE_P1_for_regression_recoded[[col]][
    Reponse_EBE_P1_for_regression_recoded[[col]] == "Do not wish to answer"
  ] <- NA
  
  # Convert to factor with levels "No" and "Yes" only (NA stays as missing)
  Reponse_EBE_P1_for_regression_recoded[[col]] <- factor(
    Reponse_EBE_P1_for_regression_recoded[[col]],
    levels = c("No", "Yes")
  )
}

Reponse_EBE_P1_for_regression_recoded$Marital_status <- 
  factor(Reponse_EBE_P1_for_regression_recoded$Marital_status,
                            levels = c("Single", "Couple", "Married or civil union", "Widowed"))


Reponse_EBE_P1_for_regression_recoded <- Reponse_EBE_P1_for_regression_recoded %>%
  mutate(Tutoring_participation_influence = factor(Tutoring_participation_influence,
                                                levels = c("Very negative", "Negative", "Neutral", "Positive", "Very positive"),
                                                ordered = F))

Reponse_EBE_P1_for_regression_recoded <- Reponse_EBE_P1_for_regression_recoded %>%
  mutate(Private_preparation_influence = factor(Private_preparation_influence,
                                                   levels = c("Very negative", "Negative", "Neutral", "Positive", "Very positive"),
                                                   ordered = F))

# Create a named vector for mapping
translation_map <- setNames(Recode_of_every_possible_answer$translation_num,
                            Recode_of_every_possible_answer$translation)

# Fonction de recodage qui garde les valeurs non trouvées
safe_recode <- function(x, map) {
  recoded <- map[as.character(x)]
  ifelse(is.na(recoded), x, recoded)
}

Reponse_EBE_P1_for_regression_recoded <- Reponse_EBE_P1_for_regression_recoded %>%
  mutate(across(
    where(~ is.character(.x) && all(!is.na(suppressWarnings(as.numeric(.x))))),
    ~ as.numeric(.x)
  ))



Reponse_EBE_P1_for_regression_recoded <- Reponse_EBE_P1_for_regression_recoded %>%
  mutate(
    MBI_SS_burnout_status = if_else(
      Emotional_Exhaustion > 14 & Cynicism > 6 & Academic_Efficacy < 23,
      1,
      0
    )
  )




univariate_lm_analysis <-
  function(data, var_to_explain, explanatory_vars){
    explanatory_vars %>%       # begin with variables of interest
      str_c(paste0(var_to_explain, " ~ "), .) %>%         # combine each variable into formula ("outcome ~ variable of interest")
      # iterate through each univariate formula
      map(                               
        .f = ~lm(                       # pass the formulas one-by-one to glm()
          formula = as.formula(.x),      # within glm(), the string formula is .x
          data = data)) %>%       # dataset
      # tidy up each of the glm regression outputs from above
      map(
        .f = ~tbl_regression(
          .x)) %>%
      tbl_stack()
  }

multivariate_lm_analysis_m1AgeGender <- function(data, var_to_explain, explanatory_vars) {
  explanatory_vars %>%
    map(function(var) {
      formula_str <- paste0(var_to_explain, " ~ Age + Gender + Current_track + ", var)
      model <- lm(formula = as.formula(formula_str), data = data)
      
      tbl_regression(
        model,
        include = all_of(var)  # show only the main tested variable
      )
    }) %>%
    tbl_stack()
}


multivariate_lm_analysis_m2_forWarwick_Edinburg_total <- function(data, var_to_explain, explanatory_vars) {
  explanatory_vars %>%
    map(function(var) {
      formula_str <- paste0(var_to_explain, " ~ Academic_Efficacy + Cynicism + Emotional_Exhaustion + SCOFF_total + ", var)
      model <- lm(formula = as.formula(formula_str), data = data)
      
      tbl_regression(
        model,
        include = all_of(var)  # show only the main tested variable
      )
    }) %>%
    tbl_stack()
}


multivariate_lm_analysis_m2_for_Academic <- function(data, var_to_explain, explanatory_vars) {
  explanatory_vars %>%
    map(function(var) {
      formula_str <- paste0(var_to_explain, " ~ Warwick_Edinburg_total + Cynicism + Emotional_Exhaustion + SCOFF_total + ", var)
      model <- lm(formula = as.formula(formula_str), data = data)
      
      tbl_regression(
        model,
        include = all_of(var)  # show only the main tested variable
      )
    }) %>%
    tbl_stack()
}
multivariate_lm_analysis_m2_forcynism <- function(data, var_to_explain, explanatory_vars) {
  explanatory_vars %>%
    map(function(var) {
      formula_str <- paste0(var_to_explain, " ~ Academic_Efficacy + Warwick_Edinburg_total + Emotional_Exhaustion + SCOFF_total + ", var)
      model <- lm(formula = as.formula(formula_str), data = data)
      
      tbl_regression(
        model,
        include = all_of(var)  # show only the main tested variable
      )
    }) %>%
    tbl_stack()
}



multivariate_lm_analysis_m2_foremotionalexhaustion <- function(data, var_to_explain, explanatory_vars) {
  explanatory_vars %>%
    map(function(var) {
      formula_str <- paste0(var_to_explain, " ~ Academic_Efficacy + Cynicism + Warwick_Edinburg_total + SCOFF_total + ", var)
      model <- lm(formula = as.formula(formula_str), data = data)
      
      tbl_regression(
        model,
        include = all_of(var)  # show only the main tested variable
      )
    }) %>%
    tbl_stack()
}



univariate_clm_analysis <-
  function(data, var_to_explain, explanatory_vars){
    explanatory_vars %>%       # begin with variables of interest
      str_c(paste0(var_to_explain, " ~ "), .) %>%         # combine each variable into formula ("outcome ~ variable of interest")
      # iterate through each univariate formula
      map(                               
        .f = ~clm(                       # pass the formulas one-by-one to glm()
          formula = as.formula(.x),      # within glm(), the string formula is .x
          data = data)) %>%       # dataset
      # tidy up each of the glm regression outputs from above
      map(
        .f = ~tbl_regression(
          .x,
          exp = T)) %>%
      tbl_stack()
  }


multivariate_clm_analysis_m1AgeGender <- function(data, var_to_explain, explanatory_vars) {
  explanatory_vars %>%
    map(function(var) {
      formula_str <- paste0(var_to_explain, " ~ Age + Gender + Current_track + ", var)
      model <- clm(formula = as.formula(formula_str), data = data)
      
      tbl_regression(
        model,
        include = all_of(var),  # inclure uniquement la variable d’intérêt
        exponentiate = TRUE
      )
    }) %>%
    tbl_stack()
}




multivariate_clm_analysis_m2Burnout_wellbeing <-
  function(data, var_to_explain, explanatory_vars) {
  explanatory_vars %>%
    map(function(var) {
      formula_str <- paste0(var_to_explain, " ~ Academic_Efficacy + Cynicism + Emotional_Exhaustion + Warwick_Edinburg_total + ", var)
      model <- clm(formula = as.formula(formula_str), data = data)
      
      tbl_regression(
        model,
        include = all_of(var),  # inclure uniquement la variable d’intérêt
        exponentiate = TRUE
      )
    }) %>%
    tbl_stack()
}

clm_univariateSCOFF <-
  univariate_clm_analysis(data = Reponse_EBE_P1_for_regression_recoded, 
                        var_to_explain = "SCOFF_total", 
                        explanatory_vars =
                          colnames(
                            Reponse_EBE_P1_for_regression_recoded %>%
                              dplyr::select(-c("SCOFF_total"))))

## Warning: (1) Hessian is numerically singular: parameters are not uniquely determined 
## In addition: Absolute convergence criterion was met, but relative criterion was not met

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

clm_SCOFF_m1 <-
  multivariate_clm_analysis_m1AgeGender(data = Reponse_EBE_P1_for_regression_recoded, 
                        var_to_explain = "SCOFF_total", 
                        explanatory_vars =
                          colnames(
                            Reponse_EBE_P1_for_regression_recoded %>%
                              dplyr::select(-c("SCOFF_total", "Age", "Gender", "Current_track"))))

## Warning: (1) Hessian is numerically singular: parameters are not uniquely determined 
## In addition: Absolute convergence criterion was met, but relative criterion was not met

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

clm_SCOFF_m2 <-
  multivariate_clm_analysis_m2Burnout_wellbeing(data = Reponse_EBE_P1_for_regression_recoded, 
                        var_to_explain = "SCOFF_total", 
                        explanatory_vars =
                          colnames(
                            Reponse_EBE_P1_for_regression_recoded %>%
                              dplyr::select(-c("SCOFF_total", "Academic_Efficacy", 
                                               "Cynicism",
                                   "Emotional_Exhaustion", "Warwick_Edinburg_total"))))

## Warning: (1) Hessian is numerically singular: parameters are not uniquely determined 
## In addition: Absolute convergence criterion was met, but relative criterion was not met

## ! `broom::tidy()` failed to tidy the model.

## ✔ `tidy_parameters()` used instead.

## ℹ Add `tidy_fun = broom.helpers::tidy_parameters` to quiet these messages.

tbl_SCOFF_merged <- 
  tbl_merge(
  tbls = list(
    clm_univariateSCOFF,
    clm_SCOFF_m1,
    clm_SCOFF_m2
  ),
  tab_spanner = c(
    "**Univariate clm (SCOFF)**",
    "**Multivariate clm m1: + Age + Gender + Current_track**",
    "**Multivariate clm m2: + Burnout + Well-being**"
  )
)


tbl_SCOFF_merged %>%
  as_tibble() %>%
  writexl::write_xlsx("./Result/tbl_SCOFF_merged.xlsx")

Reponse_EBE_P1_for_regression_recoded$Warwick_Edinburg_total <-
  as.numeric(Reponse_EBE_P1_for_regression_recoded$Warwick_Edinburg_total )

lm_univariateWE <-
  univariate_lm_analysis(data = Reponse_EBE_P1_for_regression_recoded, 
                        var_to_explain = "Warwick_Edinburg_total", 
                        explanatory_vars =
                          colnames(
                            Reponse_EBE_P1_for_regression_recoded %>%
                              dplyr::select(-c("Warwick_Edinburg_total"))))


lm_multivariateWE_m1 <-
  multivariate_lm_analysis_m1AgeGender(data = Reponse_EBE_P1_for_regression_recoded, 
                        var_to_explain = "Warwick_Edinburg_total", 
                        explanatory_vars =
                          colnames(
                            Reponse_EBE_P1_for_regression_recoded %>%
                              dplyr::select(-c("Warwick_Edinburg_total", "Age", "Gender", "Current_track"))))


lm_multivariateWE_m2 <-
  multivariate_lm_analysis_m2_forWarwick_Edinburg_total(data = Reponse_EBE_P1_for_regression_recoded, 
                        var_to_explain = "Warwick_Edinburg_total", 
                        explanatory_vars =
                          colnames(
                            Reponse_EBE_P1_for_regression_recoded %>%
                              dplyr::select(-c("SCOFF_total", "Academic_Efficacy", 
                                               "Cynicism",
                                   "Emotional_Exhaustion", "Warwick_Edinburg_total"))))

tbl_WE_merged <-
  tbl_merge(
  tbls = list(
    lm_univariateWE,
    lm_multivariateWE_m1,
    lm_multivariateWE_m2
  ),
  tab_spanner = c("**Univariate lm (Warwick Edinburg)**",
                  "**Multivariate lm m1: + Age + Gender + Current_track**",
                  "**Multivariate lm m2: + Burnout + TCA**")
)

tbl_WE_merged %>%
  as_tibble() %>%
  writexl::write_xlsx("./Result/tbl_WE_merged.xlsx")

# S'assurer que la variable est numérique
Reponse_EBE_P1_for_regression_recoded$Emotional_Exhaustion <-
  as.numeric(Reponse_EBE_P1_for_regression_recoded$Emotional_Exhaustion)

# Univarié
lm_univariateEE <-
  univariate_lm_analysis(data = Reponse_EBE_P1_for_regression_recoded, 
                         var_to_explain = "Emotional_Exhaustion", 
                         explanatory_vars =
                           colnames(
                             Reponse_EBE_P1_for_regression_recoded %>%
                               dplyr::select(-c("Emotional_Exhaustion"))))

# Multivarié m1 : ajusté pour Age + Gender + Current_track
lm_multivariateEE_m1 <-
  multivariate_lm_analysis_m1AgeGender(data = Reponse_EBE_P1_for_regression_recoded, 
                                       var_to_explain = "Emotional_Exhaustion", 
                                       explanatory_vars =
                                         colnames(
                                           Reponse_EBE_P1_for_regression_recoded %>%
                                             dplyr::select(-c("Emotional_Exhaustion", "Age", "Gender", "Current_track"))))




# Multivarié m2 : ajusté pour bien-être & autres dimensions du burnout
lm_multivariateEE_m2 <-
  multivariate_lm_analysis_m2_foremotionalexhaustion(data = Reponse_EBE_P1_for_regression_recoded, 
                                                        var_to_explain = "Emotional_Exhaustion", 
                                                        explanatory_vars =
                                                          colnames(
                                                            Reponse_EBE_P1_for_regression_recoded %>%
                                                              dplyr::select(-c("SCOFF_total", "Academic_Efficacy", 
                                                                               "Cynicism",
                                                                               "Emotional_Exhaustion", "Warwick_Edinburg_total"))))

# Tableau fusionné
tbl_EE_merged <-
  tbl_merge(
    tbls = list(
      lm_univariateEE,
      lm_multivariateEE_m1,
      lm_multivariateEE_m2
    ),
    tab_spanner = c("**Univariate lm (Emotional Exhaustion)**",
                    "**Multivariate lm m1: + Age + Gender + Current_track**",
                    "**Multivariate lm m2: + Burnout + Well-being**")
  )

tbl_EE_merged %>%
  as_tibble() %>%
  writexl::write_xlsx("./Result/tbl_EE_MBISS_merged.xlsx")

# S'assurer que la variable est numérique
Reponse_EBE_P1_for_regression_recoded$Cynicism <-
  as.numeric(Reponse_EBE_P1_for_regression_recoded$Cynicism)

# Univarié
lm_univariateCY <-
  univariate_lm_analysis(data = Reponse_EBE_P1_for_regression_recoded, 
                         var_to_explain = "Cynicism", 
                         explanatory_vars =
                           colnames(
                             Reponse_EBE_P1_for_regression_recoded %>%
                               dplyr::select(-c("Cynicism"))))

# Multivarié m1 : ajusté pour Age + Gender + Current_track
lm_multivariateCY_m1 <-
  multivariate_lm_analysis_m1AgeGender(data = Reponse_EBE_P1_for_regression_recoded, 
                                       var_to_explain = "Cynicism", 
                                       explanatory_vars =
                                         colnames(
                                           Reponse_EBE_P1_for_regression_recoded %>%
                                             dplyr::select(-c("Cynicism", "Age", "Gender", "Current_track"))))


  




# Multivarié m2 : ajusté pour bien-être & burnout
lm_multivariateCY_m2 <-
  multivariate_lm_analysis_m2_forcynism(data = Reponse_EBE_P1_for_regression_recoded, 
                                                        var_to_explain = "Cynicism", 
                                                        explanatory_vars =
                                                          colnames(
                                                            Reponse_EBE_P1_for_regression_recoded %>%
                                                              dplyr::select(-c("SCOFF_total", "Academic_Efficacy",  "Cynicism",  "Emotional_Exhaustion", "Warwick_Edinburg_total"))))

# Tableau fusionné
tbl_CY_merged <-
  tbl_merge(
    tbls = list(
      lm_univariateCY,
      lm_multivariateCY_m1,
      lm_multivariateCY_m2
    ),
    tab_spanner = c("**Univariate lm (Cynicism)**",
                    "**Multivariate lm m1: + Age + Gender + Current_track**",
                    "**Multivariate lm m2: + Burnout + Well-being**")
  )

tbl_CY_merged %>%  as_tibble() %>%
  writexl::write_xlsx("./Result/tbl_cY_MBISS_merged.xlsx")

# Transformation en numérique
Reponse_EBE_P1_for_regression_recoded$Academic_Efficacy <-
  as.numeric(Reponse_EBE_P1_for_regression_recoded$Academic_Efficacy)

# Analyse univariée
lm_univariate_AE <- univariate_lm_analysis(
  data = Reponse_EBE_P1_for_regression_recoded,
  var_to_explain = "Academic_Efficacy",
  explanatory_vars = colnames(
    Reponse_EBE_P1_for_regression_recoded %>%
      dplyr::select(-"Academic_Efficacy")
  )
)

# Multivariée m1
lm_multivariate_AE_m1 <- multivariate_lm_analysis_m1AgeGender(
  data = Reponse_EBE_P1_for_regression_recoded,
  var_to_explain = "Academic_Efficacy",
  explanatory_vars = colnames(
    Reponse_EBE_P1_for_regression_recoded %>%
      dplyr::select(-c("Academic_Efficacy", "Age", "Gender", "Current_track"))
  )
)

# Multivariée m2
lm_multivariate_AE_m2 <- multivariate_lm_analysis_m2_for_Academic(
  data = Reponse_EBE_P1_for_regression_recoded,
  var_to_explain = "Academic_Efficacy",
  explanatory_vars = colnames(
    Reponse_EBE_P1_for_regression_recoded %>%
      dplyr::select(-c("SCOFF_total", "Academic_Efficacy", 
                       "Cynicism", "Emotional_Exhaustion", "Warwick_Edinburg_total"))
  )
)

# Fusion des tableaux
tbl_AE_merged <- tbl_merge(
  tbls = list(
    lm_univariate_AE,
    lm_multivariate_AE_m1,
    lm_multivariate_AE_m2
  ),
  tab_spanner = c("**Univariate lm (Academic Efficacy from 'French MBI-SS')**",
                  "**Multivariate lm m1: + Age + Gender + Current_track**",
                  "**Multivariate lm m2: + Burnout + Wellbeing**")
)

tbl_AE_merged  %>%  as_tibble() %>%
  writexl::write_xlsx("./Result/tbl_AE_MBISS_merged.xlsx")

# EE >14 and CY > 6 and professional/academic efficacy (the equivalent of PA) <23 

# Prepare data
burnout_summary <- Reponse_EBE_P1_for_regression_recoded %>%
  count(MBI_SS_burnout_status) %>%
  mutate(
    status_label = ifelse(MBI_SS_burnout_status == 1, "Burnout", "No Burnout"),
    pct = n / sum(n) * 100
  )

# Create stacked barplot
ggplot(burnout_summary, aes(x = "", y = pct, fill = status_label)) +
  geom_bar(stat = "identity", width = 0.5) +
  geom_text(aes(label = paste0(round(pct, 1), "%")), 
            position = position_stack(vjust = 0.5), 
            color = "white", size = 5) +
  scale_fill_manual(values = c("No Burnout" = "steelblue", "Burnout" = "darkorange")) +
  labs(
    title = "Proportion of Participants by Burnout Status (MBI-SS)",
    x = NULL,
    y = "Percentage",
    fill = "Burnout Status"
  ) +
  my_theme +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank()) +
    coord_flip()

univariate_binomial_analysis <- function(data, var_to_explain, explanatory_vars) {
  explanatory_vars %>%
    str_c(paste0(var_to_explain, " ~ "), .) %>%
    map(~ glm(
      formula = as.formula(.x),
      data = data,
      family = binomial(link = "logit")
    )) %>%
    map(~ tbl_regression(
      .x,
      exponentiate = TRUE
    )) %>%
    tbl_stack()
}


multivariate_binomial_analysis_m1AgeGender <- function(data, var_to_explain, explanatory_vars) {
  explanatory_vars %>%
    map(function(var) {
      formula_str <- paste0(var_to_explain, " ~ Age + Gender + Current_track + ", var)
      model <- glm(formula = as.formula(formula_str), data = data, family = binomial(link = "logit"))
      
      tbl_regression(
        model,
        include = all_of(var),
        exponentiate = TRUE
      )
    }) %>%
    tbl_stack()
}


multivariate_binomial_analysis_m2SCOFF_wellbeing <- function(data, var_to_explain, explanatory_vars) {
  explanatory_vars %>%
    map(function(var) {
      formula_str <- paste0(var_to_explain, " ~ SCOFF_total + Warwick_Edinburg_total + ", var)
      model <- glm(formula = as.formula(formula_str), data = data, family = binomial(link = "logit"))
      
      tbl_regression(
        model,
        include = all_of(var),
        exponentiate = TRUE
      )
    }) %>%
    tbl_stack()
}


glm_univariate_burnout <- univariate_binomial_analysis(
  data = Reponse_EBE_P1_for_regression_recoded,
  var_to_explain = "MBI_SS_burnout_status",
  explanatory_vars = colnames(
    Reponse_EBE_P1_for_regression_recoded %>%
      dplyr::select(-"MBI_SS_burnout_status")
  )
)

## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

glm_multivariate_burnout_m1 <- multivariate_binomial_analysis_m1AgeGender(
  data = Reponse_EBE_P1_for_regression_recoded,
  var_to_explain = "MBI_SS_burnout_status",
  explanatory_vars = colnames(
    Reponse_EBE_P1_for_regression_recoded %>%
      dplyr::select(-c("MBI_SS_burnout_status", "Age", "Gender", "Current_track"))
  )
)

## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## suppression des ex-aequos de 'x'

## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

## Warning: glm.fit: l'algorithme n'a pas convergé

## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

## Warning: glm.fit: l'algorithme n'a pas convergé

## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

glm_multivariate_burnout_m2 <- multivariate_binomial_analysis_m2SCOFF_wellbeing(
  data = Reponse_EBE_P1_for_regression_recoded,
  var_to_explain = "MBI_SS_burnout_status",
  explanatory_vars = colnames(
    Reponse_EBE_P1_for_regression_recoded %>%
      dplyr::select(-c("MBI_SS_burnout_status", "SCOFF_total", "Warwick_Edinburg_total"))
  )
)

## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

## Warning: glm.fit: l'algorithme n'a pas convergé

## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

## Warning: glm.fit: l'algorithme n'a pas convergé

## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

## Warning: glm.fit: l'algorithme n'a pas convergé

## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

## Warning: glm.fit: l'algorithme n'a pas convergé

## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1
## Warning: glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

tbl_burnout_merged <- tbl_merge(
  tbls = list(
    glm_univariate_burnout,
    glm_multivariate_burnout_m1,
    glm_multivariate_burnout_m2
  ),
  tab_spanner = c(
    "**Univariate glm (Burnout status)**",
    "**Multivariate glm m1: + Age + Gender + Current_track**",
    "**Multivariate glm m2: + SCOFF + Wellbeing**"
  )
)

tbl_burnout_merged %>%
  as_tibble() %>%
  writexl::write_xlsx("./Result/tbl_burnout_glm_merged.xlsx")

sessionInfo()

## R version 4.3.3 (2024-02-29 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 22631)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=French_France.utf8  LC_CTYPE=French_France.utf8   
## [3] LC_MONETARY=French_France.utf8 LC_NUMERIC=C                  
## [5] LC_TIME=French_France.utf8    
## 
## time zone: Europe/Paris
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ordinal_2023.12-4.1 RColorBrewer_1.1-3  GGally_2.2.1       
##  [4] ggtext_0.1.2        scales_1.3.0        broom_1.0.8        
##  [7] UpSetR_1.4.0        ggforce_0.4.2       gtsummary_1.7.2    
## [10] lubridate_1.9.3     forcats_1.0.0       stringr_1.5.1      
## [13] dplyr_1.1.4         purrr_1.0.2         readr_2.1.5        
## [16] tibble_3.2.1        ggplot2_3.5.1       tidyverse_2.0.0    
## [19] tidyr_1.3.1         writexl_1.5.0       readxl_1.4.3       
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6         bayestestR_0.16.0    xfun_0.51           
##  [4] bslib_0.9.0          insight_1.3.0        lattice_0.22-5      
##  [7] numDeriv_2016.8-1.1  tzdb_0.4.0           vctrs_0.6.5         
## [10] tools_4.3.3          generics_0.1.3       datawizard_1.1.0    
## [13] ucminf_1.2.2         pkgconfig_2.0.3      Matrix_1.6-5        
## [16] gt_0.10.1            lifecycle_1.0.4      compiler_4.3.3      
## [19] farver_2.1.2         munsell_0.5.1        htmltools_0.5.8.1   
## [22] sass_0.4.9           yaml_2.3.10          pillar_1.10.1       
## [25] jquerylib_0.1.4      MASS_7.3-60.0.1      broom.helpers_1.15.0
## [28] cachem_1.0.8         nlme_3.1-164         commonmark_1.9.5    
## [31] ggstats_0.9.0        tidyselect_1.2.1     digest_0.6.35       
## [34] stringi_1.8.3        labeling_0.4.3       labelled_2.12.0     
## [37] polyclip_1.10-6      fastmap_1.1.1        grid_4.3.3          
## [40] colorspace_2.1-0     cli_3.6.2            magrittr_2.0.3      
## [43] utf8_1.2.4           withr_3.0.2          backports_1.5.0     
## [46] timechange_0.3.0     rmarkdown_2.29       gridExtra_2.3       
## [49] cellranger_1.1.0     hms_1.1.3            evaluate_1.0.3      
## [52] haven_2.5.4          knitr_1.50           parameters_0.26.0   
## [55] markdown_1.12        rlang_1.1.3          gridtext_0.1.5      
## [58] Rcpp_1.0.12          glue_1.7.0           tweenr_2.0.3        
## [61] xml2_1.3.6           effectsize_0.8.7     rstudioapi_0.16.0   
## [64] jsonlite_1.8.8       R6_2.6.1             plyr_1.8.9